Hidden zero-temperature bicritical point in the two-dimensional anisotropic 
Heisenberg model: Monte Carlo simulations and proper finite-size scaling 
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By considering the appropriate finite-size effect, we explain the connection between Monte Carlo 
simulations of two-dimensional anisotropic Heisenberg antiferromagnet in a field and the early renor- 
malization group calculation for the bicritical point in 2 + e dimensions. We found that the long 
length scale physics of the Monte Carlo simulations is indeed captured by the anisotropic nonlinear 
a model. Our Monte Carlo data and analysis confirm that the bicritical point in two dimensions is 
Heisenberg-like and occurs at T = 0, therefore the uncertainty in the phase diagram of this model 
is removed. 
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I. INTRODUCTION 

The two-dimensional anisotropic Heisenberg antifer- 
romagnet has recently been re-studied with extensive 
Monte Carlo simulations. 1 Twenty five years after the 
first attempt to delineate its phase diagram with Monte 
Carlo simulations^ many features of it have been clar- 
ified. However, the nature of the "spin-flop transition" 
has been an issue under debate, because the data from 
the simulations have been inadequate to show unambigu- 
ously the thermodynamic limit of the phase boundaries. 
There is no spin-flop transition in the thermodynamic 
limit according to the renormalization group (RG) cal- 
culations in 2 + e dimensions Instead, there are two 
neighboring second order phase boundaries and a disor- 
dered phase in between. By tracing the phase boundaries 
from high to low temperatures, only an upper bound for 
the bicritical temperature can be claimed. Below this 
temperature, an apparent first order spin-flop transition 
is indeed observed in both Monte Carlo simulations and 
neutron scattering experiments on quasi-two-dimensional 
Heisenberg antifcrromagnetic systems with anisotropy, 6 
It has been generally agreed that the existing data are 
consistent with the RG predictions, although some fea- 
tures near the spin-flop line, e.g. the apparent hysteresis^ 
and the unexpected crossing points in the Binder cumu- 
lant^ have not been accounted for. 

In this paper, we study the "spin-flop transition" of 
the XXZ model defined by the Hamiltonian 



SIS*) 



S*S*]-HJ2S*. (1) 



Here Sf, S, 



and Sf 



are three components of a unit vec- 
tor located on site i of a square lattice with periodic 
boundary conditions in both directions. The anisotropy 
is given by the parameter A, and H is the external mag- 
netic field in the z direction. By the term "spin-flop 
transition" , we refer to the boundary region between the 
antiferromagnetic (AF) phase and the XY phase where 




FIG. 1: Three candidates for the phase diagram. Solid lines 
are second order phase boundaries, and the dashed line in (a) 
is a first order phase boundary. 



two separate second order phase boundaries cannot be 
identified in simulations of finite-size systems. We set 
,7=1 and A = 4/5 for simplicity. Our results are ex- 
pected to be valid for < A < 1, since no qualitatively 
different behavior has been found for other values of A, 
and the phase diagram for A = 4/5 is most accurately 
known i 

At low temperatures and in a low magnetic field, sys- 
tems described by Eq. (JJJ exhibit an Ising-like AF phase, 
where the order parameter, i.e. the staggered magneti- 
zation in the z direction, has a finite value. Above a 
sufficiently large magnetic field, which is a function of 
temperature, the AF phase is replaced by an XY phase, 
where the order parameter, i.e. the in-plane staggered 
magnetization, has a finite value in the finite-size system 
due to the power-law decay of its correlation function. 
As pointed out in Ref. Ql there are three different possi- 
ble scenarios for this spin-flop transition: (a) It is a first 
order phase transition at low temperature; the first order 
phase boundary and the second order phase boundaries 
of the XY phase and the AF phase meet at a bicritical 
point, (b) The bicritical point appears at zero tempera- 
ture with a very narrow disordered phase separating the 
XY and AF phases, (c) A "biconical phase", in which 
both order parameters are nonzero, separates the XY 
and the AF phases. These three scenarios are shown in 
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Fig. ^ schematically. 

Both RG calculations-, and the Monte Carlo simula- 
tion^ have shown that (a) is realized in three dimensions, 
a scenario which is consistent with the finite critical tem- 
perature of the three dimensional isotropic Heisenberg 
model. In two dimensions, as a result of the Mermin- 
Wagner theorem^ such a bicritical point can only be at 
zero temperature, and the same prediction comes from 
RG calculations with the non-linear a modeli^^*i2*ii*i^ 
Because of the limited computer power that was then 
available, early Monte Carlo simulations with the Hamil- 
tonian of Eq. did not favor any of the above three 
scenarios^ Recently, Ref. Q] has found Ising-like scaling 
behavior in the specific heat and susceptibility on a spin- 
flop line at T/J = 0.33 for A = 4/5, which gives an 
upper bound in temperature for the bicritical point. We 
have reproduced these scaling behaviors in our simula- 
tions; however, our simulations at lower temperatures do 
not display such Ising-like scaling behavior, as we will 
show in this paper. Scenario (c) is realized in systems 
with random fields 13 - or spins with more than four com- 
ponents, e.g. the SO(5) theory>ii It is unlikely to be the 
case based on current and previous Monte Carlo simula- 
tions. 

The main issue to address in this paper is how to com- 
pare our extensive data from Monte Carlo simulations of 
limited system sizes with the RG calculations. If they are 
consistent with each other, we will have a bicritical point 
at zero temperature in the thermodynamic limit. Oth- 
erwise, we are forced to accept the apparent first order 
phase transition that we have seen at low temperatures. 
Because the correlation length is expected to diverge as 
exp(47r/T) at low temperatures, our computer power is 
very unlikely to be ample for simulations with system 
sizes larger than the correlation length in the foresee- 
able future. Therefore, properly extracting the thermo- 
dynamic limit from our data is essential. We will discuss 
the Monte Carlo simulation and derive the theoretical 
finite-size scaling predictions in Sec. ^ present results 
from our simulation and the data analysis in Sec. II I II 
and conclusions with discussions in Sec. IIVI 



continuity. In scenario (a) in Fig.l, the spin-flop line 
starts from H c at T = 0, shows a small positive slope as 
the temperature increases, and eventually develops into 
an umbilicus, where the AF-paramagnctic phase bound- 
ary and the XY-paramagnetic phase boundary are clearly 
separated. For the other possibilities in 2 dimensions, two 
phase boundaries might well be indistinguishable from 
each other but yield the same "effective" finite tempera- 
ture behavior. 

We perform Monte Carlo simulations that scan the 
magnetic field across the spin-flop transition at con- 
stant temperatures. The random spin configurations 
are generated with the heatbath algorithm,— which has 
been recently shown to be the optimal single spin algo- 
rithm for Heisenberg spins It is a rejection- free al- 
gorithm due to the fact that the probability distribution 
P(S) cx exp(h • S/T) of a single spin subject to a mag- 
netic field can be produced with uniform random number 
generators without rejecting trial flips. As in Ref. UtI 
we simultaneously perform two simulations with differ- 
ent initial configurations, the equilibrium is achieved by 
letting these two simulations run until their order pa- 
rameters are almost equal. Each simulation then runs 
for 10 6 to 10 7 Monte Carlo steps per spin at the given 
temperature and magnetic field, and is allowed to accu- 
mulate substantially more data near the critical magnetic 
field. The data for energy, magnetization, and staggered 
magnetization are stored for histogram reweighing. The 
difference in observables measured from independent sim- 
ulations are used to estimate the error bar. Totally we 
have spent about 1 CPU year and accumulated 27GB 
data, which represents a modest computing load on a 
decent Linux cluster. 

We calculate the ensemble average of the staggered 
magnetization: 



(2) 



where the sign is different on two sublattices. We define 
the tilt angle between the z axis and Mt as 



II. THEORETICAL BACKGROUND 

A. Monte Carlo simulation and traditional 
finite-size scaling 

At zero temperature, the system has a spin-flop tran- 
sition at H c — AJy/l — A 2 . The system is in the fully 
aligned AF configuration for H < H c and is in the 
spin- flop configuration for H > H c . As H — » if+, 



the spin-flop configuration approaches S* — J and 
gx, y _ ^_iy n <n,y I ^A- j where (n x ,n y ) is an arbitrary 

unit vector. Obviously, the spin-flop transition is a first 
order phase transition at zero temperature, since the or- 
der parameter (staggered magnetization) exhibits a dis- 



= cos 



Aft ' 



(3) 



where the value of 8 is restricted to [0,7r/2] due to the 
inversion symmetry of the staggered magnetization, and 
calculate its probability distribution. The Binder cumu- 
lant for an arbitrary observable X is defined as 



U*(X) = 1 



3 (A 2 )' 



(4) 



where (• • • ) denotes the ensemble average. The suscepti- 
bility per spin of Af J is defined as 



T 



{[Mtf)-{\Mt\Y 



(5) 
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and the specific heat per spin is defined as 



L 2 



(0) 



where E is the internal energy per site. 

In order to determine the type of phase transitions, we 
performed the traditional finite-size scaling analysis. In 
the case of an Ising-like second order phase transition, 
the specific heat exponent a = 0, because c v is logarith- 
mically divergent near the critical temperature or critical 
field, and its maximum value increases as In L. The local- 
ization length exponent v — 1 ; the spontaneous magneti- 
zation exponent j3 = 1/8; and the susceptibility exponent 
7 = 7/4. In case of a first order phase transition, these 
exponents can also be defined properly, and they have a 
different set of values^ a = 1,0 = 0, v = 1/2, 7=1. (In 
general v = 1/d, where d is the dimensionality of the sys- 
tem.) Near the phase transition point, it is convenient to 
define the reduced temperature t — T/T c — 1 and reduced 
field h = H/H c — 1. Finite size scaling theory predicts 
the following scaling relations: 



m?) L 

(x) L 



L^l v C(tL 1 l v ). 



(7) 
(8) 



Similar relations hold for (Mj) 2 + (M^) 2 and c v , as well 
as for transitions with h as the control variable. In par- 
ticular, for c v of Ising-like transition, the prefactor L a / V 
is replaced by C + InL, where C is a constant. 



B. Low temperature finite-size scaling predictions 

As we shall see from the results in the next section, 
neither the scaling relations for the first order phase tran- 
sition nor those for second order phase transitions fit our 
data perfectly well. The first order phase transition ap- 
parently seems to be more consistent though. To resolve 
these discrepancies, we return to the RG calculation of bi- 
critical phenomenon in 2 + e dimensions. The long length 
scale physics is expected to be governed by the non- linear 
a model with an anisotropy termi&a 



(9) 



where the 0(3) spin field (tti(x), ^(x), a) satisfies the 
constraint 7r 2 + a 2 = 1 and the anisotropy constant 
g oc H — H c . [Note: we follow the convention in 
Ref. here, so the Boltzmann factor is exp(7i cr ).] In 
general, g is expected to be a linear combination of 
both magnetic field and temperature contributions: 19 
g = a{H 2 - H 2 ) + b(T - T c ). The spin-flop line (de- 
fined by g = 0) of our model is almost horizontal, with a 
very small negative b. In the scenario shown in Fig.^b), 
H c actually falls in the paramagnetic phase between two 
phase boundaries. The RG calculation can proceed with 
either the momentum shell renormalization technique* or 



Polyakov's approach. 3,10 Same results can also be derived 
with Brezin and Zinn-Justin's methodiii*^ The RG flow 
differential equations for renormalized temperature T(l) 
and anisotropy constant g(l) are 



dT{l) 
dl 

dg(l) 
dl 



= -eT(0 
= 2 5 (0- 



n — 3 



, nif 

TQ) g(l) 
7T 1+3(0' 



1 



1 + 9(1) 



(10) 
(11) 



where e = d — 2. For our present purpose, it is sufficient 
to set e = 0, n = 3 and keep the lowest order terms of g 
in these equations. The approximate solution then reads 



T(l) = 



2ttT 



(12) 



2-k-TV 

g(l) = ge 2l T 2 /T(lf = ge 2l [l-Tl/(2*)] 2 . (13) 

Additionally, the equation for the renormalization con- 
stants for 7r and a can be derived from Ref. 0: 



din Cot 
dl 

dlnCr 
dl 



= d- 
= d- 



T(l)g(l) 

47T 5 (0 
2vr ' 



1' 



(14) 
(15) 



Using Eq. I|12|) and d = 2, for negligible g, one finds that 
Ctt and Co- are both given by 



Tl 
2^ 



(16) 



These results imply that for g > 0, the renormalization 
flow goes to g — +oo, which corresponds to the XY 
model; for g < 0, it goes to g — — oo, which corresponds 
to the two-dimensional Ising model; for g — 0, it goes 
to the high temperature paramagnetic phase, which is 
consistent with the Mermin- Wagner theorem £ Conven- 
tionally, the thermodynamic quantities are calculated by 
integrating the renormalization flow equation to a char- 
acteristic scale I*, for which g(l*) = 0(1). The cor- 
responding thermodynamic quantities can be calculated 
by other means, e.g. perturbation theory. Then a tra- 
jectory integral "matching" formalismi^fliSiiSa is used to 
obtain these quantities at the original scale I — 0. For ex- 
ample, the susceptibility has the scaling form %(T, g) ~ 
l T <&(ge^ l T ) for small T and g in the thermodynamic 
limit^ which is analogous to the crossover phenomenon of 
a finite temperature bicritical pointiSiSi The Ising or XY 
critical behavior is contained in the scaling function <E>, 
which is calculated from the effective Ising or XY model 
(depending on the sign of 3) at the scale I* ~ e 4Tr / T . Due 
to the finite size of the systems in our simulations, the 
above RG flow cannot go very far in the RG flow diagram 
for an infinite system, but only to the scale set by the sys- 
tem size V ln(L/a), where a is the microscopic lattice 
constant (fixed to be unity). For our model, L = 100a 
is not a large system size compared to the correlation 
length, so we expect the lowest order solutions Eqs. i|12|) 
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and (|13fl are sufficient. In the real space RG terminol- 
ogy^ at scale I', the spins in a block of linear size e l are 
effectively replaced by a single block spin, therefore, we 
expect our simulation with system size L in a magnetic 
field sufficiently close to the critical field is governed by 
the simple Hamiltonian 



To calculate the specific heat, the standard approach 
is to derive it from the free energy / = — L~ 2 \nZ. As 
in Ref. 0, we perform a trajectory integral from I = to 
1 = 1* to evaluate the free energy. For the clarity of our 
discussion, we reproduce the trajectory integral formula 
from Ref. J here: 



K 



eff 



-qa 2 , 



where the coefficient q is given by 



ffQnL) 
2T(lnL) 



2T 



1 - 



TlnL 



2tt 



(17) 



(18) 



Near the spin-flop transition, g = k(T,L)h + 0(h 2 ), 
where h is the reduced field, and we expect the coef- 
ficient k(T,L) to weakly depend on T and L. As a 
result, the free energy should appear as a function of 
hT- x L 2 [\ - T\nL/(2n)} 3 . One can define a "correla- 
tion length exponent" v — 1/2, based on the leading L 2 
dependence. The staggered magnetization also has to 
contain a factor from the spin renormalization constant: 

({ M lf) L = L ~ 4 C ( ff2 )w eff ' (19) 
({Ml y ) 2 ) L = {{Mlf + {Ml) 2 ) L = L-\l <7r 2 >„£2Q) 



where (. . . ) w denotes the thermal average with respect 
to the effective Hamiltonian Eq. (fP7|) . When Cr/Cr ~ 1; 
both of them have the scaling form: 



TlnL 
2tt 



•Fi 



11,-1 



hL 2 ( TlnL 

— v-— 



(21) 



Furthermore, the distribution of should be on a 
sphere of radius 1 — TlnL/(27r); the distribution of its 
tilt angle has the form 



P(0) = Y~ x exp(-q cos 2 9) sin 9, 



(22) 



where Y is a normalization constant. At the critical 
field g = h = 0, the above distribution is uniform. 
Consequently the critical Binder cumulant for cos 9 is 
Ui(cos9) = 2/5. Ui{Ml) should also be close to this 
value, with a small negative correction due to the longi- 
tudinal fluctuation of M^. The maximum susceptibility 
corresponds to the maximum fluctuation in c, so the fol- 
lowing scaling formula is expected: 



( 2 

S<7 



X 



hL 2 
T 



1 - 



TlnL 
2tt 



(23) 



One can again define a susceptibility exponent 7=1 
based on the leading L 2 dependence, which is assumed 
to be L 1 I V . The critical exponent ratio is then identical 
to that for first order phase transition; therefore, the log- 
arithmic corrections are required to determine whether 
the spin-flop transition is a first order phase transition 
or not. 



f(T,g) 



l G (l)dl 



-dl* 



f[T(l*),g(l*)}, (24) 



in which we set d = 2 in our calculation and the kernel 
Go depends on both the on-shell Green's function and 
the spin renormalization factors. With the final scale I* = 
InL and the approximate solution Eqs. I|12l) and l|13l) . this 
trajectory integral can be evaluated analytically with the 
technique in Appendix A of Ref. l22t Among many terms 
in the result, we are particularly interested in those of 
the form L _2 / S [T(Z*), g(l*)], for reasons which will soon 
become clear. It turns out that the only such term is 
fs[g(l*)} = - Ml + g(l*)]. The last term in Eq. (JUJl is 
given by the single spin Hamiltonian Eq. (|17fl : 



L~ 2 F 



gl? ( TlnL 
~2T~ V 2~^~ 



where 



F(q) = - In / exp(-g cos 2 9)2tt sin 9d9. 
Jo 

The divergent term in the specific heat is given by 
d 



dT 



T 2 L -2d(F + fs) 



dT 



(25) 



(26) 



(27) 



One would expect the leading divergent term in Eq. Ij27(l 
is proportional to L 2 g 2 , with multiplicative logarithmic 
corrections. However, obviously on the spin-flop line 
g = 0, this term vanishes. This dilemma is solved by not- 
ing that T in Eq. l|27Jl is the real temperature at which 
the simulation is performed, while T in Eqs. (|24|l and l|25|) 
is a renormalized temperature for the long length scale 
effective Hamiltonian, and from now on it will be denoted 
T* . The spin-flop line does not follow a constant temper- 
ature, in other words, the effective anisotropy g is also 
a function of T. If we change T in the simulation, we 
actually change the effective anisotropy g, which could 
drive the system to cross the spin-flop line. Therefore, 
the partial derivative in Eq. 1)2 7|l should be written as 







d__ 

df ~ a df~* 



dg 



(28) 



Two derivatives d/dg acting on L~ 2 F would produce 
a divergent term proportional to L 2 [l — T* lnL/(27r)] 6 . 
Similarly, f s results in a divergent term proportional to 
L 2 [l - T* lnL/(27r)] 4 . After all these considerations, the 
divergent part of the specific heat is expected to have the 
following scaling form: 



c v = L 2 x 6 C 6 (hL 2 x 3 ) + L 2 x 4 C 4 {hL 



(29) 
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where x = 1 — T* lnL/(27r). If the spin-flop line is per- 
fectly horizontal, i.e. 6 = 0, there will not be a peak in 
the specific heat. We expect the C§ term to be dominant 
in Eq. I|29fl . because F also contributes a factor T~ 2 after 
being differentiated twice, on the other hand, a straight- 
forward calculation shows that the C4 term does not show 
a peak at g = 0, but an asymmetric background. 

The above trajectory integral actually connects the in- 
termediate renormalized Hamiltonian Eq. © to the final 
renormalized Hamiltonian Eq. (|17(l . In doing so, we have 
actually ignored a similar trajectory integral which con- 
nects the bare Hamiltonian Eq. (JIJ to the intermediate 
Hamiltonian Eq. (JJJJ. Although we cannot write down an 
analytic expression for its integrand and its final contri- 
bution to the free energy, we expect it to depend on the 
lattice constant, which is the smallest length scale in the 
system, and an intermediate length scale at which Eq. 
is valid, but not on the largest length scale L. Therefore, 
it cannot possibly give rise to an L 2 divergence. For the 
divergent terms in c v , it is hence justified to ignore this 
precursory trajectory integral. 

Unlike the conventional finite-size scaling, the above 
finite-size scaling relations are not expected to hold for 
very large system sizes, because for sufficiently large L, 
the phase boundaries of the AF and XY phases become 
discernible. This obviously happens when TlnL ~ 2n, 
i.e. the spin renormalization constants become very 
small, corresponding to a disordered phase. For either 
large I or large g, the higher order terms in Eqs. IjlOjl and 
become non-negligible; therefore, the approximate 
solution Eqs. I|12fl and i|13fl and the above finite-size scal- 
ing relations are no longer sufficient. Thus, we expect the 
above finite-size scaling analysis are valid for small effec- 
tive anisotropy, i.e. H k H c , and system sizes smaller 
than the correlation length. 



III. DATA AND ANALYSIS 

A. Raw data and the apparent first order phase 
transition 

We first present some representative raw data for the 
staggered magnetization from simulations with different 
system sizes, and we will then show traditional finite-size 
scaling plots. Figure [21 shows the staggered magnetiza- 
tion at T = 0.2 for magnetic fields on both side of the 
finite lattice spin-flop transition. The spin-flop transition 
can be clearly observed, and when viewed at moderate 
resolution (see the insets) the data suggest a lst-order 
transition. Both ((Afj) 2 ) and ((Mj^) 2 ) reach either zero 
or a large value quickly as H deviates from H c , which is 
why we have to zoom in on a very narrow range (0.005 
in width) of the magnetic field to see the finite-size ef- 
fect for relatively large systems. The critical magnetic 
field that separates the XY phase from the AF phase can 
already be estimated rather accurately from the size de- 
pendence seen in Fig. [21 There is no noticeable sign of 
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FIG. 2: (color online) Staggered magnetization across the fi- 
nite lattice spin-flop line at T = 0.2. Each group of nearby 
data points are from one simulation performed with single 
fixed magnetic field. Histogram reweighing was used to cal- 
culate observables in slightly different magnetic fields. The 
insets show the data measured from a larger range of mag- 
netic field. 



an intermediate phase which has zero expectation values 
for both order parameters. There is also no sign of a re- 
gion with nonzero values for both order parameters. In 
the XY phase, we see from the inset of Fig. |2fb) that 
the saturated values of the order parameter are smaller 
for larger systems, which is consistent with the property 
of the XY phase. (We expect the saturated value to 
decrease as L~ 2r >, where 77 is a temperature dependent 
cxponenti^SiSl) On the other side, in the AF phase, the 
inset of Fig. [2fa) shows that ((Mj) 2 ) calculated with 
different system sizes quickly converges to the nonzero 
thermodynamic limit. 

If the spin-flop transition contains an Ising-like tran- 
sition, then we would expect to see the universal value 
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FIG. 3: (color online) Binder cumulant for various sizes at 
T = 0.2, the crossing points of these curves are near H — 
2.4055, U& at the crossing point is different from the universal 
value for Ising-like transition. 



of the Binder cumulant XJ\ ~ 0.618 at the critical mag- 
netic field, therefore we plot Ua{MI) in Fig. [3] to see if 
it indicates an Ising-like second order phase transition. 
However, the curves in Fig. [21 cross each other at values 
close to 0.4, clearly below the Ising universal value. There 
is not any systematic trend indicating that the crossing 
point moves up with increasing system size. This result 
is in agreement with Ref. Q, where it is shown that only 
when T > 0.4 do the crossing points in the Binder cu- 
mulant curves start to move towards the Ising universal 
value. As Ref.Qhas shown in the phase diagram, the AF- 
paramagnetic boundary and XY-paramagnetic boundary 
are clearly separated for T > 0.4. We have also calcu- 
lated the Binder cumulant for T — 0.1,0.265, and 0.33, 
and seen results similar to Fig. |3J So it is indeed tempt- 
ing to consider the possibility of a bicritical point which 
exists between T = 0.33 and T = 0.4. 

To see whether our data really fit the behavior ex- 
pected for a first order phase transition, we show the 
finite-size scaling analysis of M\ in Fig. According to 
Eq. Q, all the data points should collapse onto a single 
curve if we choose the right exponents v and (3. The field 
dependence of the data in both Fig. [21and Fig.|21indicate 
that H c is between 2.405 and 2.406. In Fig. g| we have 
fine tuned H c to collapse the data as much as possible 
onto a single curve, and found that it is impossible to do 
so with Ising exponents. With first order exponents, the 
scaling plot Fig 01 a) is clearly better. However, system- 
atic deviations are discernible in the low field AF phase. 
In fact, if we allow (3 to have a small value, we have found 
that the scaling plot looks best with (3 = 0.031. The same 
phenomenon has also been observed in traditional finite- 
size scaling plots for ((AfJ y ) 2 ). 

Similar scaling analyses have been done for T = 0.1, 
and 0.265. In all those cases, the Ising exponents fail, but 
the first order exponents work reasonably well, especially 
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FIG. 4: (color online) Finite-size scaling of Mj. (a) Expo- 
nents used belong to first order phase transition; (b) Ising 
exponents are used. Error bars are not plotted since most of 
them are smaller or equal to the size of the symbols, as shown 
in Fig. |3 



for T — 0.1. If we accept a small which goes to zero as 
T goes to zero, then at finite temperatures, the sponta- 
neous staggered magnetization does decay as M| oc \h\P 
in the thermodynamic limit (L — > oo). We would indeed 
have a second order phase transition. However the corre- 
lation length exponent v — 1/2 is not consistent with this 
scenario. For susceptibility and specific heat, finite-size 
scaling analysis has also ruled out Ising-like transition as 
a possible scenario. Especially for specific heat, its max- 
imum value scales roughly as L 2 at T = 0.1, 0.2, which is 
glaringly different from the Ising-like logarithmic behav- 
ior. 



B. Proper finite-size effect and scaling 
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FIG. 5: Probability distribution of the staggered magnetiza- 
tion across the spin-flop transition. P(r, z) is proportional to 
the grey scale. The data are from simulations with L = 100 
at T = 0.265. 

To characterize the spin-flop transition, we first study 
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FIG. 6: (color online) Probability distribution of the tilt angle 
9. Symbols are calculated from the histogram of 9 in the 
simulations, and the solid lines are curve fitting with Eq. (1221 . 
The fitting parameter q is plotted in the inset, where the 
straight line is a linear fit. The simulations were performed 
with L = 80 at T = 0.265. 



FIG. 7: (color online) Finite-size scaling analysis correspond- 
ing to Eq. at T = 0.2. (a) The AF order parameter; 
(b) the XY order parameter. We introduce T* as a fitting 
parameter which is slightly above the actual temperature at 
which the simulation is done. Error bars are smaller or equal 
to the size of the symbols. 



the probability distribution of the staggered magnetiza- 
tion to get a vivid picture of the spin-flop transition. 
Due to the obvious rotational symmetry around the z 
axis, we plot in Fig.Elthe probability distribution P(r, z) 
of systems with L = 100 in cylindrical coordinates at 
T = 0.265. It can be clearly seen from this figure that 
is distributed on a sphere of about 0.65 in radius. In 
the AF phase, the distribution is confined to the north 
and south poles; in the XY phase, the distribution is 
near the equator. One would naturally define a finite- 
size critical field H c at which the distribution is uniform 
on this sphere, i.e. zero effective anisotropy. For L = 100 
and T — 0.265, we have found H c « 2.411 which is in 
agreement with the observation of ((MJ) 2 ) and Ui(M%). 
At this critical field, J7 4 (M]) « 2/5, which can be imme- 
diately calculated given this uniform distribution. The 
same pattern of transition have been observed at different 
temperatures with different system sizes. For larger sys- 
tems and higher temperatures, we have observed smaller 
spheres. Clearly, a finite-size effect can be extracted from 
these spherical distributions. 

To be more quantitative, we calculate the angular dis- 
tribution of the staggered magnetization. Figure shows 
the distribution of the tilt angle 9 calculated from simu- 
lations at T = 0.265 for L = 80. The data are fitted with 
Eq. I|22l) . which only has one fitting parameter q. Obvi- 
ously this simple functional form fits the simulations very 
well. The inset in Fig.|B|shows the fitting parameter, i.e. 
Eq. i|18|) . is a linear function of H, which is consistent 
with our assumption that g oc h. Thus, we have shown 
that close to H c , the thermodynamics of the staggered 
magnetization is indeed governed by the simple Hamil- 
tonian Eq. l(T7|> . 

The staggered magnetization clearly has a finite value 



at H c in the simulations. If this is also true for the ther- 
modynamic limit, then we will have to choose between 
the scenario of a first order phase transition, or that of 
an intermediate ordered phase with a tetracritical point. 
However, this is not the case for two-dimensional sys- 
tems. Suppose we enlarge the simulation to a size nL, 
since we have shown that an L x L block indeed be- 
haves as an effective Heisenberg spin in an anisotropic 
potential, the long length scale physics should be cap- 
tured by a Hamiltonian of these block spins. Although 
microscopically, the anisotropy may have different forms, 
at a large enough length scale the go 1 term turns out 
to be the only dominant term, as shown by the proba- 
bility distribution in Fig. In terms of RG terminol- 
ogy, it is the only relevant perturbation that keeps the z 
axis rotational symmetry. Since a similar anisotropy in 
the kinetic energy, gi(d fJl a) 2 has been found to be irrele- 
vant^ our simulations have indeed justified Eq. @ as an 
appropriate Hamiltonian to analyze the spin-flop transi- 
tion. Following the RG analysis outlined in Sec. [HI one 
would conclude that the behavior in the thermodynamic 
limit is two second order phase boundaries and a bicrit- 
ical point at zero temperature, i.e. that which is shown 
in Fig. [Tib). 

Now we perform a finite-size scaling analysis in order 
to examine whether or not the logarithmic corrections in 
Eq. (|21|l can be seen in the data from the simulations. 
The renormalized temperature T in the effective Hamil- 
tonian Eq. @ is, in general, different from the tempera- 
ture at which the simulation is performed. Therefore, we 
are allowed to use an extra fitting parameter T* in Fig.[7| 
to make the data points collapse better on a single curve. 
We find that all the data points do collapse onto a single 
curve in both Fig. Ufa) and[7Jb) with the choice T* = 0.3. 
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Actually, if this renormalized temperature is not intro- 
duced, Fig. [7| still looks clearly better than Fig. Ufa), 
where we assume a first order phase transition. In the 
transition region, we also expect ((A/]) 2 + (Mj ) 2 ) to be 
nearly a constant as seen in Fig. [5] This is also consistent 
with Fig. where the decrease in ((Mj) 2 ) is compen- 
sated by the increase in ^(Mj y ) 2 ). A small /? improves 
Fig. EI a) because a logarithmic correction enters by the 
Taylor expansion for small (3: L 2 ®l v « 1 + 2/3z/~ 1 mZ. 
Error bars are not plotted in Fig. and the other scaling 
plots, because most of them are smaller or equal to the 
size of the symbols used in these figures. They can be 
found in Fig. [2] and |3J Figure [S] shows the same scaling 
analysis for simulations performed at T = 0.1. We have 
found that this scaling plot is very clear even without in- 
troducing a renormalized temperature. The critical field 
H c is fixed up to six significant digits in Fig. [7| and [H] by 
closely examining the central critical region of very small 
h. However, one should not push this too far because 
the apparent existence of a single critical field is only a 
finite-size effect. For large systems, there are two sepa- 
rate critical fields for the AF phase boundary and the XY 
phase boundary respectively. The same finite-size scaling 
analysis at T = 0.265 still produces a very good single 
curve; however, the best scaling is achieved by allowing 
H c for the XY phase to be slightly above that for the AF 
phase, behavior which could be understood as a sign of 
separating phase boundaries. 

The logarithmic correction for 1VF comes from the spin 
renormalization constants in Eqs. I|19|l and i|20|) . The 
effective spin a and 7T can be calculated directly from 
M'f, given that Q ~ ( a . One would expect that the 
finite-size scaling of (<r 2 ) only involves a scaling of the 
magnetic field. In fact, Eqs. (|17l) and (|18|l imply that 
(er 2 ) L for different sizes as well as different temperatures 
can collapse on to a single curve, provided that fc(T, L) — 
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FIG. 8: (color online) The same as Fig. [T] but for simulations 
at T = 0.1. Unlike Fig.|7| we have used T* =T here. 
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FIG. 9: (color online) Finite-size scaling plot of the effective 
a spin component, (a) from simulations at T = 0.2, (b) from 
simulations at T — 0.1. Instead of adjusting T* , we have 
simply set T* = T. (a) and (b) are expected to displays the 
same curve, see text for details. Error bars are smaller or 
equal to the size of the symbols. 
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FIG. 10: (color online) Scaling plot for susceptibility near the 
spin-flop transition according to Eq. 1231 . (a) at T = 0.1, (b) 
at T = 0.2. 



h/g depends weakly on T and L. Figure plots (c 2 ) L 
versus hT~ l Li 2 {l — TlnL/27r) 3 , and shows that this is 
indeed true in the simulations. We have also calculated 
the Binder cumulant ^(cos^). Since the longitudinal 
fluctuation in is completely projected out, U^{cos9) 
is exactly 2/5 at the crossing point, which corresponds 
to the critical magnetic field for finite-size systems. 

As for the susceptibility, we test Eq. (|2*3t by plot- 
ting xL~ 2 (l - TlnL/27r)- 2 versus hL 2 (l - TlnL/2ir) 3 
in Fig. El For L — 40, 60, and 80, our data collapse 
onto single curves at T = 0.1 and T = 0.2 very well. For 
two larger sizes, L = 100 and L = 120, although we have 
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FIG. 11: (color online) Scaling plot for specific heat near the 
spin-flop transition according to Eq. 12911 . (a) at T = 0.1, (b) 
at T = 0.2. 

fewer data points near the peak of the susceptibility and a 
few data points have larger statistical errors than others, 
they appear to fall on the same curve reasonably well. We 
have also plotted the data assuming finite-size scaling for 
an Ising-like transition or a first order phase transition 
but found none of those works better than Fig. I1UI 

Finally, in Fig.^2 we show the finite-size scaling plots 
for specific heat which follow Eq. 129fl . Only the first Cq 
term has been considered here. Most of the data points 
collapse quite well onto a single curve in both Fig. Illf al 
and II If , and this verifies the validity of our finite-size 
scaling analysis. We have made the same scaling plot for 
specific heat at T = 0.265 and 0.33, and found none of 
them obeys this scaling formula. This is obviously be- 
cause at higher temperatures, phase boundaries of the 
XY and AF phases start to separate. In fact, those data 
at higher temperatures are consistent with an Ising tran- 
sition. As in other figures, a group of nearby data points 
are calculated from the same simulation with histogram 
reweighing. The disadvantage of histogram reweighing is 
that the statistical error in one simulation is replicated by 
several data points. Therefore, the deviations observed 
in Fig. El as weu as Fig. are n °t real "systematic" 
deviations, but an artifact of histogram reweighing. 

IV. DISCUSSION 

Although near the critical magnetic field, the simula- 
tion is hindered by the huge correlation length, so that 
two separate second order phase boundaries may never 
be revealed, the staggered magnetization behaves as a 
rigid spin of nearly fixed length subject to an anisotropic 
potential. The length of the spin, and the form and 
strength of the anisotropy are predictions of the RG cal- 
culation, which have been confirmed by our Monte Carlo 
simulations. Logarithmic corrections of particular forms, 



which are absent in either first order phase transition 
or second order phase transitions, have been found in 
the Monte Carlo simulations in complete agreement with 
our theoretical predictions. For the phase diagram of the 
two-dimensional XXZ Heisenberg antiferromagnet with 
an easy axis, we have thus reached the conclusion that 
in the thermodynamic limit, the bicritical point appears 
at T = 0, with an exponentially narrow paramagnetic 
phase sandwiched between the low field AF phase and 
the high field XY phase. The location of the bicritical 
point is "hidden" from detection by "ordinary" means 
because of finite-size effects that must be carefully an- 
alyzed. In the simulation of any finite-size system, the 
staggered magnetization remains non-zero on the "spin- 
flop line" which only becomes two separate phase bound- 
aries in the thermodynamic limit. The Ising-like or XY- 
like critical behavior near the spin-flop transition is also 
expected to be seen in an exponentially narrow range of 
magnetic field. Our results should be valid for all two- 
dimensional Heisenberg antiferromagnet or ferromagnet 
with short range interaction and easy axis anisotropy of 
different forms, since their long length scale physics is 
described by the same Hamiltonian. Due to the finite 
experimental resolution, including effects from the inho- 
mogeneity in magnetic field and disorder in samples, and 
possible crossover to three dimensional behaviors, we ex- 
pect experiments would only reveal apparent first order 
spin-flop transitions and an apparent bicritical point at 
finite temperature. Also symmetry breaking perturba- 
tions, e.g. crystal fields of square symmetry, are highly 
relevant to the bicritical phenomena in 2+e dimensions as 
well as the intricate correlations in the XY phaseA2& As 
a result, both the AF phase and the "XY" phase might 
have discrete symmetries, so that the two second order 
phase boundaries are likely to be reduced to a first order 
phase boundary. 



V. CONCLUSION 

We have carried out extensive Monte Carlo simula- 
tions for a two-dimensional anisotropic Heisenberg anti- 
ferromagnet near the spin-flop transition and performed 
a finite-size scaling analysis based on the anisotropic non- 
linear a model in the regime where the correlation length 
is much larger than the system size. We have found that, 
although the finite-size effects tend to obscure the asymp- 
totic behavior, the results of the Monte Carlo simula- 
tions are perfectly consistent with the RG calculation for 
the spin-flop transition and that the bicritical point is at 
T = 0. 
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